Pairing in graphene: A quantum Monte Carlo study 
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To address the issue of electron correlation driven superconductivity in graphene, we perform 
a systematic quantum Monte Carlo study of the pairing correlation in the t — U — V Hubbard 
model on a honeycomb lattice. For V = and close to half filling, we find that pairing with 
d + id {d^2_y2 + id'^y in its specific form) symmetry dominates pairing with extended-s symmetry. 
However, as the system size or the on-site Coulomb interaction increases, the long-range part of the 
d + id pairing correlation decreases and tends to vanish in the thermodynamic limit. An inclusion of 
nearest-neighbor interaction V , either repulsive or attractive, has a small effect on the extended-s 
pairing correlation, but strongly suppresses the d + id pairing correlation. 



Recently, graphene has attracted the attention of ex- 
perimentahsts and theoristsi"— . One of the most in- 
triguing properties of graphene is that its chemical po- 
tential can be tuned through an electric field effect, 
and hence it is possible to change the type of carriers, 
electrons, or holes, opening the doors for carbon based 
electronica^i^. Doped graphene has a finite density of 
state at the chemical potential, which, in combination 
with pronounced antiferromagnetic (AFM) spin fluctu- 
ations close to half fillingli^, may lead to an unconven- 
tional superconductivity. Experimentally, superconduct- 
ing (SC) states in graphene have been realized by the 
proximity effect through contact with SC electrodes^, 
which indicates that Cooper pairs can propagate coher- 
ently in graphene. These facts raise the question as to 
whether it would be possible to modify graphene to be 
an intrinsic superconductor. 

Various theoretical attempts^"— have been made to 
understand the superconductivity in graphene. Uchoa et 
alJ^ suggested that an extended-s (ES) SC phase may be 
realized at the mean-field level due to the special struc- 
ture of the honeycomb lattice. On the other hand, in a 
weak-coupling functional renormalization group studj^, 
Honerkamp found that with a nearest-neighbor (NN) 
spin-spin interaction J, doping away from half filling can 
lead to a. d + id SC state, which is similar to the SC 
state on the triangular latticei^iii. This d + id SC state 
was also found to be stable in a mean-field studyi^ of 
a phenomenological Hamiltoniar*i^. Recent variational 
Monte Carlo simulations of the repulsive Hubbard model 
provide further support for the d + id SC statei^. 

Although the results based on mean-field theory and 
other approximate methods are encouraging, it is far 
from certain that there exists a SC ground state in the 
physical parameter region of graphene. It is well known 
that the low energy properties of graphene can be de- 
scribed by the two-dimensional Hubbard model on a hon- 
eycomb lattice^. In graphene, the on-site Hubbard repul- 
sion is approximately half the band width, and this places 
graphene in an intermediate-coupling regime. Thus, it is 
questionable to approach the effect of electron correla- 



tions in graphene from either a weak-coupling or strong- 
coupling limit, as was done in many previous theoretical 
studies. In view of the above-mentioned facts, we em- 
ploy two accurate numerical methods, i.e., the determi- 
nant quantum Monte Carlo (DQMC)^ and constrained 
path Monte Carlo (CPMC) methodsi^i^a, to investigate 
possible electron correlation driven superconductivity in 
graphene. 

Our extensive numerical simulations show that close 
to half filling, pairing with d + id symmetry is dominant 
over pairing with ES symmetry. However, the long-range 
part of the d + id pairing correlation tends to vanish in 
the thermodynamic limit, suggesting the absence of elec- 
tron correlation driven superconductivity in our studied 
model. We also find that the NN interaction, either re- 
pulsive or attractive, does not enhance the tendency to 
the d + id or ES superconductivity. 

The structure of graphene can be described in terms 
of two interpenetrating triangular sublattices, A and B, 
and its low-energy electronic and magnetic properties can 
be well described by the extended Hubbard model on a 
honeycomb lattico^, 

H = -t^al„bi+,,c^ +h.c. + U'^{naifnaii + nbi'fnbii) 

iqa i 

+ V^nainM+jj + fi^iriaia +nbia), (1) 
iri ia 

Here, (c^L) annihilates (creates) electrons at site 
with spin a {cr=^, 1) on sublatticc A, (^Icr) annihilates 
(creates) electrons at the site Ri with spin a {cr—'[, 1) on 
sublattice B, Uaia = and nMc = b\^ha- t w 2.5eV 

is the NN hopping integral and /i the chemical potential. 
U w 6eV and V denote the on-site Hubbard interaction 
and NN interaction, respectively. We have mainly used 
U = 3|t| and = in this Rapid Communication, except 
as explicitly noted otherwise. 

Our numerical calculations are performed on lattices of 
double-48, double-75, double-108, and double-147 sites 
with periodic boundary conditions. The double-48 lat- 
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FIG. 1: (Color online) (a) Sketch of graphene with double-48 
sites, (b) Phases of the d + id and ES pairing symmetries 
on the honeycomb lattice. For clarity, the phase for the ES 
symmetry is omitted, which equals for all Si. 



tice is sketched in Fig. [TJa), where blue circles and yel- 
low circles indicate the A and B sublattices, respectively. 
The system is simulated using DQMC at finite temper- 
ature and CPMC at zero temperature. The basic strat- 
egy of DQMC is to express the partition function as a 
high-dimensional integral over a set of random auxiliary 
fields. The integral is then accomplished by Monte Carlo 
techniques. In the CPMC method, the ground-state 
wave function is projected from an initial wave function 
by a branching random walk in an ovcrcomplctc space 
of constrained Slater determinants, which have positive 
overlaps with a known trial wave function. Extensive 
benchmark calculations showed that the systematic er- 
ror induced by constraint is within a few percent and the 
ground-state observablcs are insensitive to the choice of 
trial wave functioni^. In our CPMC simulations, we em- 
ploy closed-shell electron fillings and use the correspond- 
ing free-electron wave function as the trial wave function. 

As magnetic excitation might play an important role 
in the SC mechanism of electronic correlated systems, we 
first study the magnetic correlations in graphene. Specif- 
ically, we compute the NN spin correlation jy = 
{Sf ■ Sj) and the spin structure factor S'(q), which is 
defined as. 



(mM-mj,,), (2) 



d.d'—a^h ij 



where m^^ = al^Ui^ - al^Uii and nii^ = bl^bi^-bl^bii. Ns 
denotes the number of sublattice sites. 

To investigate the SC property of graphene, we com- 
pute the pairing susceptibility. 



and the pairing correlation, 

C„(r = Ri-Rj) = (At«AJj)), 



(3) 



(4) 



where a stands for the pairing symmetry. Due to the 
constraint of the on-site Hubbard interaction in Eq. (1), 
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FIG. 2: (Color online) (a) Spin structure factor S{q) on the 
double-48 lattice along the high symmetry lines in the first 
Brillouin zone (BZ) at electron fillings < n 1.00, 1.10, 
and 1.20. (b) NN spin correlation S^i^^ versus < n > on 
different lattices. The inset of (a) shows the first BZ, with a 
denoting the distance between NN lattice sites. The results 
are presented for temperature r=|t|/6. 



pairing between two sublattices is favored and the corre- 
sponding order parameter Ajj(i) is defined as 



(5) 



with being the form factor of pairing function. 

Here, the vectors 5i {I = 1, 2, 3) denote the NN intersub- 
lattice connections, as sketched in Fig.[IJb). Considering 
that the pairing symmetry of graphene is governed by 
the D6 point group, two form factors for NN pairing de- 
scribed by the Ai and E2 irreducible representations of 
the D6 point group arc given bjsii. 
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In Fig. [21 we present the spin structure factor S{q) and 
the NN spin correlation at different electron fill- 

ings for temperature r=|t|/6. A broad peak between M 
and K in Fig. [H^a) and a negative NN spin correlation in 
Fig. [2Ib) indicate the existence of AFM spin correlation 
in graphene close to half filling. From Fig. [H^a), one can 
notice that even at half-filling, the peak is rather weak, 
suggesting that magnetic ordering is strongly frustrated 
due to the structure of the honeycomb lattice. Moreover, 
Fig. [5] shows that as the electron filling < n > increases 
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FIG. 3: (Color online) (a) Pairing susceptibility Pa as a func- 
tion of temperature T for different pairing symmetries and 
different fillings, (b) Pd+id and Pd+id as a function of tem- 
perature at < n >= 1.05 and < n >= 1.10. Dashed and 
dashed-dotted lines in (a) represent the d + id pairing suscep- 
tibility at [/ = 0. The inset of (b) shows the effective pairing 
interaction Pd+id — Pd+id as a function of temperature. 
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FIG. 4; (Color online) Pairing correlation Ca as a function 
of distance for the d + id and ES pairing symmetries on the 
double-75 lattice with < n >= 1.107 (a) and double-108 lat- 
tice with < n >= 1.13 (b). The inset of (b) shows the average 
of the long-range d + id pairing correlation Cd+id vs -^j^ for 
different f/'s at a filling < n >~ 1.1. 



from half filling, ^(q) is reduced in the region around 
the K point and j> becomes less negative, which 
indicates that the AFM spin correlation is suppressed 
when the system is doped away from half filling. As it 
is expected that fermion systems with strong on-site re- 
pulsion may exhibit superconductivity induced by AFM 
spin fluctuations, and from the behavior of magnetic cor- 
relation shown in Fig. [21 it seems that the electron cor- 
relation driven superconductivity is possible through a 
similar mechanism in graphene. In the following, we dis- 
cuss the behavior of pairing susceptibility and pairing 
correlation in the low doping region. 

Figure [3] shows the temperature dependence of pairing 
susceptibilities for different pairing symmetries and elec- 
tron fillings on the double-48 lattice. From Fig. ^a.), it 
is clear to see that within the filling range investigated, 
the pairing susceptibilities for both d+id and £'5' pairing 
symmetries increase as the temperature is lowered. Most 
remarkable is that Pd+id increases much faster than Pes 
at low temperatures. This demonstrates that the d + id 
pairing symmetry is dominant over the ES pairing sym- 
metry in the low doping region. In the whole tempera- 
ture regime, one can observe that the value of Pd+id at 
U = 3\t\ is smaller than the corresponding noninteract- 
ing one {U = 0), as displayed in Fig. [3l^a), which reflects 
the fact that the reduction of quasiparticlc weight (self- 



energy effect) due to electron correlation plays a negative 
role in enhancing the pairing susceptibility. 

In order to extract the effective pairing interaction 
in different pairing channels, the bubble contribution 
Pa{i,j) is also evaluated, which is achieved by replac- 
ing (a,Vjta!+54^j+5,/t> with (a^^&jt)(al+5,i^J+'5,'t) in 
Eq. (3). In Fig. El (b), we plot both Pd+id and Pd+^d 
for comparison. It is apparent that Pd+id shows a very 
similar temperature dependence to that of Pd+id- The 
effective pairing interaction, which can be estimated by 
the difference between Pd+id and Pd+id, is found to take 
a positive value and to increase with lowering tempera- 
ture, as clearly shown in the inset figure. The positive 
effective pairing interaction indicates that there actually 
exists attraction for the d + id pairing. 

Based on the DQMC results for the pairing susceptibil- 
ity, one expects that there may exist the d + id SC state 
in the low-tcmpcrature region, which is manifested by a 
divergence of Pd+id at a certain temperature. Unfortu- 
nately, it is not clear whether the pairing susceptibility 
keeps growing at low temperatures since the sign prob- 
lem prevents simulation in the low-temperature regime. 
In order to shed light on the critical issue as to whether 
there exists a long-range off-diagonal d + id SC order in 
the ground state, we now turn to discuss the results ob- 
tained from the CPMC method. 
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FIG. 5: (Color online) Pairing correlation as a function of 
distance for the d+id (filled symbols) and ES (open symbols) 
pairing symmetries on the double-75 lattice with < n >= 
1.107. The value of NN interaction V is indicated by the shape 
of the symbol, which is also applied for the ES symmetry. 



In Fig. |4l we compare the long-range part of pairing 
correlations for the d + id and ES pairing symmetries on 
the double- 75 and double- 108 lattices, with a filling of 
approximately < n >= I.l. One can readily see that 
Cd+idir) is larger than CesIi^) for all long-range dis- 
tances between electron pairs. Similar behavior is also 
observed on the double- 147 lattice with < n >= 1.095 
(not shown here). This re-enforces our finding that the 
d + id pairing symmetry dominates the ES pairing sym- 
metry in the low doping region. 

To gain insight into the behavior of the d + id pair- 
ing correlation in the thermodynamic limit, we examine 
the evolution of Cd+id with increasing lattice size. In the 
inset of Fig. IH^b), the average of the long-range d + id 
pairing correlation, Cd+id = Sr>4a Cd+rd{r), where 
N' is the number of electron pairs with r > 4a, is plotted 
as a function of for U = 0,2|t|,3|i|, and 4|f|. We 

observe that Cd+id decreases as the lattice size increases, 
and shows a clear tendency to vanish in the thermody- 
namic limit (-^^ -^> 0) at [/ = 0. This result, together 

with a decrease of Cd+id with increasing J7, suggests the 
absence of long-range d + id SC order in the parameter 
regime investigated. Our finding is in agreement with the 
functional renormalization group studjJ^, where the SC 
instability does not occur in the doped Hubbard model. 
A similar decrease of C d+id with increasing lattice size 



was also observed in the variational Monte Carlo calcula- 
tionsi^. Thus, although AFM fluctuations can mediate a 
NN singlet pair, quantum fluctuations beyond the mean- 
fleld level shall destroy the phase coherence between elec- 
tron pairs. 

In graphene, the long-range interaction may also play 
a significant role on its physical properties, especially in 
the low doping region, where the long-range Coulomb 
interaction is not effectively screened because of a small 
density of state at the Fermi energy. We have studied the 
effect of NN interaction on the pairing correlation. In Fig. 
[5l the pairing correlations for both d + id and ES pair- 
ing symmetries are displayed as a function of distance on 
the double-75 lattice with different NN interaction V' s. 
Here, we consider both repulsive and attractive NN in- 
teractions. We notice that the ES pairing correlation is 
hardly affected by the NN interaction, whereas the d + id 
pairing correlation is suppressed by either a repulsive or 
attractive NN interaction. Our results are contrary to 
the finding of Uchoa et alt^, where the NN attraction 
can stabilize the ES SC state at the mean-field level. In 
addition, the inclusion of NN interaction does not en- 
hance the tendency to the d + id SC state. Therefore, we 
can conclude that within the extended Hubbard model, 
there seems to be an absence of SC order in the ground 
state. 

In summary, we have studied the behavior of pair- 
ing correlation within the extended Hubbard model on 
a honeycomb lattice by using quantum Monte Carlo sim- 
ulations. The results obtained from both DQMC and 
CPMC show that close to half filling, pairing with d + id 
symmetry dominates pairing with ES symmetry, which is 
consistent with previous mean-field and functional renor- 
malization group studies. This provides strong evidence 
that pairing with different spatial phases is favored for 
the AFM fluctuation mediated pairing interaction, which 
is similar to the case on the triangular lattic o^^'^^ . How- 
ever, the tendency of the d + id pairing correlation to 
vanish in the thermodynamic limit suggests that electron 
correlation in graphene is not strong enough to produce 
an intrinsic superconductivity. In an induced SC state 
by the proximity effect through a connection to another 
superconductor, the dominant d + id pairing could be 
manifested in the unique properties of the SC gap func- 
tion and the Andrecv conductance spectr 
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